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Abstract 

This numerical study examines the importance of self-consistently accounting for transport and elec- 
^ ■ trostatics in the calculaiton of semiconductor/metal Schottky contact resistivity. It is shown that ignoring 
^ ■ such self-consistency results in significant under-estimation of the contact resistivity. An explicit numerical 
pLn '• method has also been proposed to efficiently improve contact resistivity calculations. 
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I. Introduction 

In modern MOSFET designs, silicon/silicide Schottky hetero-interfaces are commonly used to form 
¥ ; source/drain contacts. As the device size continues to shrink with each generation, the resistance associated 
?h ■ with the silicon/silicide contacts begins to have its impact on the device performance [1]. Therefore, it 
g . becomes important to achieve an accurate yet efficient method to evaluate silicon/silicide contact resistivity 
\ at different doping concentrations and temperatures. Some widely refered earlier work on J-V relation 
; and contact resistance of Schottky contact are due to [2], [3], [4]. They identified three regimes where the 
Q major current contribution comes from field emission, thermionic field emission or thermionic emission, 
^ ! depending on doping concentration and temperature. Analytical expressions of contact resistivity were 
£j | obtained for each of the three regimes. However, their works were based on two assumptions. Firstly, 
Poisson equation was only solved by assuming fully depletion of free carriers within the depletion region, 
i.e. the band-bending of silicon near the interface is parabolic. Secondly, the continuity equation was not 
! solved, and therefore the carrier quasi-Fermi level was assumed to be flat across the depletion region. To 
>■ ; examine the validity of those two assumptions, we have implemented in our device simulator Prophet a 
■ self-consistent Schottky barrier diode simulation model proposed by [5]. In this physical model, Poisson 
and continuity equations are simultaneously solved with the distributed tunneling current self-consistently 
included. By comparing analytical models against this physical model, we are able to demonstrate that 
CN ; those two assumptions lead to significant deviation in contact resistivity calculations, particularly at high 
qq ■ doping concentration of modern device and high temperature under ESD condition. In this work, we also 
O ■ extend the analytical model by removing those two assumptions. In the improved model, a unified contact 
j> \ resistance expression is explicitly obtained^, and the results show excellent match with those from the 
physical model. 
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II. Models and Results 

In Fig.l is a schematic energy band-diagram of a silicon/silicide Schottky contact. In this work, we 
assume the majority carriers are electrons. The hetero-interface in this plot is at location w and the 
depletion region is between and w. The difference of silicon/silicide affinity is q<p b where q is the 
elementary charge. The doping density in silicon region is N D and the effective density of states is N c . 
The difference between conduction band and quasi-fermi level at charge neutral region is q<p s . The applied 
forward bias is Vf. If we only consider the contact resistance at very low bias, we have Vf — > 0. The 
barrier height E b can be expressed as E b = q((p b — <fr s — Vf). We base this work on Maxwell-Boltzmann 
statistics for mathematical simplicity H. All the derivations below can be readily extended to Fermi-Dirac 



'Numerical integrations are involved in the expression. 

2 Also because the physical model is currently implemented based on M-B statistics. 
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statistics, and the effects investigated in this work should still be present in that case. It can be shown 
that exp(—q(j) s /kT) = N D /N C under M-B statistics. 

Following the treatment of [6], [7], the general tunneling probability at loaction x within the depletion 
region is given by 



t{E(x)) = exp[ — —J ^2m*[E(x') - E(x)]dx' ^ 

= exp ( / y/2m*[E(x') - E(x)]——rdx') (1) 
v n J e(x) aa,\x ) j 



' E(x) 

where m* is the tunneling effective mass, E(x) is the conduction bandedge energy at position x. The 
electro- static potential (f>(x) is given by 4>(x) = —E(x)/q. It can be seen that, if dx/dE(x) is known, 
the tunneling probability r(E) can be evaluated explicitly by performing the numerical integral in Eqn. 
1. Since the bandedge energy E(x) is a monotonic function of position x, it is convenient to use a 
dimensionless quantity a = E(x)/E b as the basic variable [3]. In doing so, the general expression for the 
tunneling probability is re-written as 



where E$q = (qh/ATc)^jNr>/m*€. The function F(a) carries the information of potential variation in the 
depletion region and is defined as 

f< ) - e ( dE \ 2 - eEb ( da ) 2 m 

{a> ~ 2q 2 N D E b [ dx ' 2q 2 N D { dx ' ' K ) 

In the literatures, different energy band-diagrams (E(x)) have been assumed to calculate r(a). In the 
work of [3], they based their calculation completely on the parabolic band-bending relation 

E(x) = qN D x 2 /{2e) (4) 

and obtained an analytical expression 

r c - R -{a)=exp[-^y{a)l (5) 

where y(a) = \/l — a — alog[(l + \/l — ot)/y/a\. However, the work in [2], [5] was based on another 
widely used formular for tunneling probability: 



t p - s -(E) = exp 



8tt V2~^F (E b - Ef/ 2 
1, h IFI 



(6) 



where F is the electrical field perpendicular to the interface. Matsuzawa et al. used |F(x)| = (E b 
E(x))/(w — x) in their work [5], and therefore the tunneling probability can be re- written as 



r RS -(E) = exp 



8nV2 



E b -E-(w-x) . (7) 
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This expression can also be directly derived from the general expression (Eqn. [B by assuming linear 
relation between E and x. It should be noted that although triangular potential barrier was assumed in 
obtaining Eqn. [71 the position x appears explicitly in the formula, which still contains the band-bending 
information. In our physical model, the treatment follows that of [5], and Eqn. [7] is adopted. In the 
simulation, the realistic potential variation enters through the relation of x(E) by solving the Poisson 
and continuity equations. In order to have fair comparisons between numerical models and the physical 
model, we use Eqn. [7] for the numerical models in most cases throughout this work. 
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The general expression for the forward current density (electrons inject from silicon to silicide) is given 
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(8) 



where A = A*T 2 for Richardson constant A* = Anm* qk 2 / h 3 , r\ is the local electron quasi-Fermi level 
variation compared with that at charge neutral region. Therefore, if constant quasi-Fermi level is assumed, 
77 is always zero within the entire silicon region. The first term in the bracket of Eqn. [8] corresponds 
to field emission and thermionic field emission, while the second term is from the thermionic emission. 
Similarly, the reverse current density (electron inject from silicide to silicon) is given by 



Jr= kf ( 



-q<t>s/kT 



r{E)e 



E + qVf 

^ dE + e" 



(9) 



The total net current density as a function of Vf is therefore obtained as 

J = Jj — J r 
A 



-q<t> 3 /kT 



kT 



r(E)e kT (e 



<m 

kT 



~)dE + e 



IVb 
kT 



kT ) 
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where rjh is the value of rj at the interface. The contact resistivity at zero bias is then obtained as R = 
(dJ/dVf^lvf^o [4]. In previous work [2], [3], continuity equation was not solved. Therefore, the variation 
of rj within the depletion region was ignored. Under such an assumption, the total current density is 
obtained as 



J 



A r -g<t>s/kT 

kT 



r(E)e~ w r dE + e 



1-e 



qVj_ 

kT 



(ID 



Hence, the inverse of contact resistivity at zero bias is 



R- 



Aq 

Jkfy' 

AqE b 



-q<t>slkT 



-q<t>s/kT 



r{E)e"kTdE + e" 



r(a)e kT a da + e 



kT 



{kTr y o ._ ■ (12) 

In Fig. 2 is the contact resistivity obtained from simulations of the physical model for N D = le20cm 3 
and T = 300, 500, 700, 900K, respectively. We firstly compare it with results from a numerical model, 
namely model A. In model A, the parabolic potential variation is assumed by substituting Eqn. @] into 
Eqn. |7l which gives 

4 E b 



r v« ) = exp[-^— -VT 

o -£^00 



a. 



(1-v^)]. 



(13) 



In this model, quasi-Fermi level is regarded as flat in silicon, i.e. Eqn. [T2] is used. It can be seen in 
Fig. 2 that severe discrepancy of calculated resistivity exists between the numerical model A and the 
physical model throughout the entire temperature range. In order to investigate its cause, we compare the 
potential variation computed from these two models in Fig. 3. It is clearly observed that, the assumption 
of fully depletion becomes invalid at energies near the quasi-Fermi level. Since the tunneling probability 
exponentially depends on the tunneling distance, this discrepancy in x — E relation leads to significant 
difference in the tunneling probability, as shown in Fig. 4. By removing the parabolic band-bending 
assumption in model A, the accuracy of contact calculation can be greatly improved. For this purpose, 
the one-dimensional Poisson equation needs to be solved for the depletion region. For M-B statistics, the 
Poisson equation is expressed as 

ld 2 E 



q 



q dx 2 



qN D 



D 



n) 



1-e 



kT 



(14) 
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Multiply dE/dx on both sides of Eqn. [14] and integrate from to E, and we obtain 



dE 
dx 

Express it in terms of a and we have 



2(?N D f E H _jl 



1 - e~^r)dE 



o 



1/2 



(15) 



F(a) = / [l-e-^ a )da, (16) 



according to the definition of F(a) in Eqn. [3] Eqn. [16] therefore defines the potential variation from the 
exact solution to Poisson equationj- Apply Eqn. [3] and [16] to Eqn. [7J and it is obtained that 



r B (a) = exp 



2 Eb . u /2 f 1 da 
[I -a) i- 



[a H 



(17) 



3 E 00 J a y/F(t 

We then have an improved numerical model, namely model B, in which Eqn. [16] and [FT] are used to 
compute the tunneling probability, and constant quasi-Fermi level (Eqn. [Til) is still assumed. We also 
plot the x — E relation and the tunneling probability of model B in Fig. 3 and Fig. 4, respectively. They 
show excellent match with those of the physical model. The computed resistivity from model B is also 
plotted in Fig. 2. It can be seen that qualitative improvement is obtained over model A with reference 
to the physical model. At moderate temperature, the match between model B and the physical model is 
fairly good. However, an evident discrepancy between them is still observed at high temperature, which is 
addressed in the next paragraph. As we previously mentioned, our numerical model A, B and the physical 
model are all based on Eqn. [7J which partially assumed triangular potential barrier a priori. Therefore, 
it would be interesting to see what the effect of the potential variation is based on the the more general 
expression of tunneling probability (Eqn. [2]). As mentioned earlier, the expression used in [3] (Eqn. 0) 
is derived from Eqn. [2] by assuming parabolic band-bending completely. We use their expression Eqn. [5] 
in numerical model A. In another numerical model B we use the general expression Eqn. [2] and the 
realistic band-bending formula Eqn. [16] The contact resistivity computed using these two models are 
plotted in Fig. 5, and similar trend is observed: the assumption of the parabolic band-bending severely 
under-estimates the resistivity at all temperature range. 

As shown in Fig. 2, significant discrepancy in calculated resistivity exists at high temperature between 
numerical model B and the physical model. However, it can be seen from Fig. 4 that the tunneling 
probability matches for the two models. Therefore, the source of this discrepancy originates from the 
assumption of constant rj in the derivation from Eqn. [10] to Eqn. [TTJ In Fig. 6 is the band-diagram of 
a Schottky contact simulated by the physical model for forward bias Vf = 0.1V. The variation of the 
electron quasi-Fermi level within the depletion region is evident due to finite carrier supply rate by drift- 
diffusion. If we let Vf — > 0, we expect this variation i] to approach zero at the same time. But the ratio 
rj/Vf can be finite in this limit. Since the resistivity is a differential quantity, it can be affected by this 
effect. However, it should be noted that the validity of drift-diffusion model in the depletion region is an 
open question itself, since the depletion width is comparable to electron mean free path. If the carrier 
supply is not limited by the drift-diffusion process, the variation of quasi-Fermi level can be negligible. 
A possible way to examine this problem is Monte Carlo simulation. The electron continuity equation is 

^- + U = 0, (18) 

dx 



3 After assuming constant quasi-Fermi level and ignoring hole concentrations. 
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where j is electron flux due to carrier transport in the depletion region, and U is tunneling flux density. 
Plug in the tunneling probability and integrate from to to and we obtain that 
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(19) 
(20) 
(21) 



If we assume the electron transport in the depletion region can be modeled as drift-diffusion, as in the 
physical model, the electron flux can also be expressed as 

dr] 



j(a) = /J,N D e~ 



kT ' £ kT 



dx 



/jN D e k T e ^ , 



ehb da 

where fi is electron mobility. Compare Eqn. [21] and [221 and let Vf, r\ — > 0, then we obtain 

dlogil-v/Vf) _ 



(22) 



or equivalently, 



da 



rj/Vf = exp 



A _ q<t>s 

e kT 



kTq\i 



2M 



Pa) 



d e kT 



(23) 



A _ q<t>s 

e kT 



eE h 
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The inverse of resistivity near zero bias is then revised as 



! AqE± ft 
{kT) 2 



r a e kT 



[1 - rj/V f )da + e 



The resistivity calculated using this improved model (model C) is also plotted in Fig. 2, and 
match with the results of the physical model is observed. 
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Fig. 1. A schematic plot of energy diagram of a silicon/silicide Schottky contact. 




Fig. 2. 



Contact resistivity vs. temperature calculated by various numerical models and the physical model. 
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Fig. 3. Distance to the interface vs. barrier energy (x-E relation) calculated by various numerical models and the physical model. 
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Fig. 4. 



Tunneling probability vs. barrier energy calculated by various numerical models based on Eqn. [7] and the physical model. 
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Fig. 5. Contact resistivity calculated by two numerical models based on Eqn. [2] 




Fig. 6. 



Energy band diagram of a Schottky contact at 0.1V forward bias simulated by the physical model. 



